--- title: Basics keywords: fastai sidebar: home_sidebar nb_path: "nbs/01_basics.ipynb" ---
Fastai 2 provides support to read and display medical images using pydicom and pillow, however only 2D images can be read. Faimed3d provides tool to load an preocess 3d medical data
Fastai already provides PILDicom and TensorDicom. Since these objects are designed to handle 2D data, using then might lead to problems in later tasks. So a custom Tensor class will be created to be able to use some of the very handy fastai functions.
freqhist_bins and hist_scaled are needed to for image normalization
The retain_type function from fast ai can be used, to retain the TensorDicom3D type.
retain_type(t, typ = TensorDicom3D)
3D data is stored in a variety of data formats (DICOM, NIfTI, NRRD, Analyze, ...) which should be supported by the image loader. Also DICOM data is often stored as individual slices (DICOM series) and not as a volume. SimpleITK is a powerfull library which can handle many data formats, including all of the above mentioned and many more. It will therefore be used to read the 3D volumes.
im = readdcm_3d('/media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM')
im.shape, show_image_3d(im, axis = 0, nrow = 10)
Sometimes multiple 3D images (e.g. a batch) need to be displayed. With a wrapper for show_image_3d this is conveniently possible.
im.shape
show_images_3d(torch.stack((im,im)), axis = 0, nrow = 31, figsize = (25, 15))
Often the important areas of the images are also located centrally. For example in CT scans, the skin should always be included in the image, which means that the outer pixels will display only air or the scanner table. These areas can be cropped, effectivly reducing the size of the image without reducing the resolution.
cropped_im_abs = crop_3d(im, ((0,0),(100,100),(100,100)), perc_margins=False)
cropped_im_perc = crop_3d(im, ((0.,0.),(0.1,0.1),(0.1,0.1)), perc_margins=True)
im.shape, cropped_im_abs.shape, cropped_im_perc.shape
show_image_3d(cropped_im_abs)
Medical data has different resolutions. Most CT scans will be at 512x512 px but for MRI the resolution varies. Also the number of slices is different for every scan. To train a neural network all input needs to be the same size, thus a function to resize the images alongside each axis is needed.
# works well but is overly complicated. Using F.interpolate is much simpler and is implemented as a tfms.
def resize_3d_tensor(t: Tensor, new_shape: int):
"""
A function to resize a 3D image using torch.nn.functional.grid_sample
Args:
t (Tensor): a Rank 3 Tensor to be resized
new_dim (int): a tuple with the new x,y,z dimensions of the tensor after resize
Returns:
a resized `Tensor`. Note that the header with metadata is lost in this process.
Procedure:
[fake RGB] -> [batch dim] -> [flow field] -> [resampling] -> [squeeze]
fake RGB:
Create a fake RGB 3D image through generating fake color channels.
btach dim:
add a 5th batch dimension
flow field:
create a flow-field for rescaling:
a. create a 1D tensor giving a linear progression from -1 to 1
b. creat a mesh-grid (the flow field) from x,y,z tensors from (a)
Taken form the offical Pytroch documention:
Given an input and a flow-field grid, computes the output using input values and pixel locations from grid.
In the spatial (4-D) case, for input with shape (N,C,Hin,Win) and with grid in shape (N, Hout, Wout, 2), the output will have shape (N, C, Hout,Wout)
In the case of 5D inputs, grid[n, d, h, w] specifies the x, y, z pixel locations for interpolating output[n, :, d, h, w].
mode argument specifies nearest or bilinear interpolation method to sample the input pixels.
resampling:
resample the input tensor according to the flow field
squeeze:
remove fake color channels and batch dim, returning only the 3D tensor
"""
if type(new_shape) is tuple and len(new_shape) == 3:
z,x,y = new_shape # for a reason, I do currently not understand, order of the axis changes from resampling. flipping the order of x,y,z is the current workaround
else:
raise ValueError('"new_dim" must be a tuple with length 3, specifying the new (x,y,z) dimensions of the 3D tensor')
t = t.unsqueeze(0) # create fake color channel
t = t.unsqueeze(0).float() # create batch dim
x = torch.linspace(-1, 1, x) # create resampling 'directions' for pixels in each axis
y = torch.linspace(-1, 1, y)
z = torch.linspace(-1, 1, z)
meshx, meshy, meshz = torch.meshgrid((x, y, z)) #
grid = torch.stack((meshy, meshx , meshz), 3) # create flow field. x and y need to be switched as otherwise the images are rotated.
grid = grid.unsqueeze(0) # add batch dim
out = F.grid_sample(t, grid, align_corners=True, mode = 'bilinear') # rescale the 5D tensor
return out.squeeze().permute(2,0,1).contiguous() # remove fake color channels and batch dim, reorder the image (the Z axis has moved to the back...)
resized_im = resize_3d_tensor(im, (100,50,50))
resized_im.shape, show_image_3d(resized_im)
Especially for MRI the pixel values can vary between scanner types. However, the imagenet stats, as provided by fastai should probably not be used, as they are specfic for color images and not MRI images. The optimal solution would probably be to calculate the stats on the present dataset.
Statistics of intensity values, pooled across the whole training dataset (raw, unscaled pixel values and pixel values scaled between 0 and 1)
| Statistic | ADC | T2 | T1 map |
|---|---|---|---|
| Pooled Maximum of intensity values | 3064.2500 | 1369.5652 | 4095 |
| Pooled Minimum of intensity values | 0. | 0. | 0. |
| Pooled Mean of intensity values | 511.1060 | 259.1454 | 740.8268 |
| Pooled SD of intensity values | 488.8707 | 190.2448 | 688.8238 |
| Scaled Maximum of intensity values | 1. | 1. | 1. |
| Scaled Minimum of intensity values | 0. | 0. | 0. |
| Scaled Mean of intensity values | 0.1675 | 0.1918 | 0.1809 |
| Scaled SD of intensity values | 0.1599 | 0.1409 | 0.1682 |
However, just scaling wiht one mean and std might not be the optimal solution (see the excelent Kaggle kernel form Jeremy Howard why). Although for MRI images it might be ok, as the pixel values are closer together.
im1 = resize_3d_tensor(readdcm_3d('/media/ScaleOut/prostata/data/dcm/A0042280702/T2/DICOM'), (10, 100, 100))
im2 = resize_3d_tensor(readdcm_3d('/media/ScaleOut/prostata/data/dcm/A0041912240/T2/DICOM'), (10, 100, 100))
im1.shape, im2.shape
im = TensorDicom3D.create('/media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM')
im.metadata['table']
Segmentation Masks need an own class, because image transforms can only be applied to prior specified classes. Thus, transforms altering the pixel value (e.g. add noise, change brightness) will not be applied to Masks if the class differs from the original image.
class TensorMask3D(TensorDicom3D):
"Base class for 3d Segmentation, inherits from TensorDicom3D"
@classmethod
def create(cls, fn:(Path,str,Tensor,ndarray), **kwargs)->None:
"open the mask, keep as float"
if isinstance(fn,ndarray): return cls(fn)
if isinstance(fn, Tensor): return cls(fn)
# nifti files differ in pixel orientation from DICOM.
# This might be due to differnt pixel orientations.
# currently this problem was solved otherwise, so the create function does not differ from TensorDicom3D.
# However future flipping and Rotation might be implemented here.
return cls(load_image_3d(fn))
px1 = TensorDicom3D(im1.flatten())
plt.hist(px1.numpy(), bins = 40)
plt.hist(px1.hist_scaled().numpy(), bins = 40)
plt.title("Distribution of pixel values in T2 images original and scaled")
plt.show()
The pixels are now all nearly equally distributed.
im1 = TensorDicom3D(im1)
im2 = TensorDicom3D(im2)
show_images_3d(torch.stack((im1, im2, im1.hist_scaled(), im2.hist_scaled()), axis = 0), nrow = 10)
I belive both steps should executed after reading the images, even before cropping and resizing, thus they will be directly integrated into the readdcm_3d function.
normalize(im).mean(), normalize(im).std()
im = readdcm_3d('/media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM', return_scaled=True, return_normalized=True)
im.shape, show_image_3d(im, axis = 0, nrow = 10)
Somtimes the mask is better viewed as 3D object. Rendering is implemented as described in this example: https://scikit-image.org/docs/dev/auto_examples/edges/plot_marching_cubes.html A faster, more felxible way might be using ipyvolume or vtk.
To implement 3D rendering for the mask, the TensorMask3D class needs to be expanded
im = TensorMask3D.create('/media/ScaleOut/prostata/data/dcm/A0042197734/T2/Annotation/Annotation.nii.gz')
im.render_3d(alpha = (0.15, 0.15))
im.calc_volume()